PINN with Data
Solid Mechanics Example
By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH
Table of Contents
1. Thin Plate (Solid Mechanics)Ā¶
1.1. Problem SetupĀ¶
- We will solve thin plate equations to find displacement and stress distribution of thin plate
- Based on Kirchhoff-Love plate theory, three hypotheses were used
- straight lines normal to the mid-surface remain straight after deformation
- straight lines normal to the mid-surface remain normal to the mid-surface after deformation
- the thickness of the plate does not change during a deformation
- A non-uniform stretching force is applied to square elastic plate
- Only one quarter of the plate is considered since the geometry and in-plane forces are symmetric (yellow domain)
- Problem properties
$$
E = 50 \operatorname{Mpa}, \quad \nu = 0.3, \quad \omega = 20 \operatorname{mm}, \quad h = 1 \operatorname{mm}, \quad f = 1 \operatorname{Mpa}
$$
- Governing equations (Fƶpplāvon KĆ”rmĆ”n equations) for the isotropic elastic plate:
$$
\begin{align*}
&{E \over 1 - \nu^2}\left({\partial^2 u \over \partial x^2} + {1 - \nu \over 2}{\partial^2 u \over \partial y^2} + {1 + \nu \over 2}{\partial^2 v \over \partial x \partial y} \right) = 0\\\\
&{E \over 1 - \nu^2}\left({\partial^2 v \over \partial y^2} + {1 - \nu \over 2}{\partial^2 v \over \partial x^2} + {1 + \nu \over 2}{\partial^2 x \over \partial x \partial y} \right) = 0
\end{align*}
$$
- Two Dirichlet boundary conditions at $x = 0,\, y = 0\; (B.C.ā , B.C.ā”)$:
$$
v(x,y) = 0 \qquad \text{at} \quad y = 0\\\\
u(x,y) = 0 \qquad \text{at} \quad x = 0
$$
- Two free boundary conditions at $y = \omega / 2\; (B.C.ā¢)$:
$$
\sigma_{yy} = 0,\quad \sigma_{yx} = 0
$$
- Free boundary condition and in-plane force boundary condition at $x = \omega / 2\; (B.C.ā£)$:
$$
\sigma_{xx} = P \centerdot h,\quad \sigma_{xy} = 0
$$
- Make a neural network and loss funcitons like below:
1.2. Numerical SolutionĀ¶
- Numerical solution of this problem is illustrated in below figures
- $x, y$ direction displacement and stress $u$, $v$, $\sigma_x$, $\sigma_y$, respectively
- Solve this problem using PINN and then compare with a numerical solution
InĀ [1]:
import numpy as np
import matplotlib.pyplot as plt
import random
import deepxde as dde
import pandas as pd
import torch
from sklearn.metrics import r2_score
InĀ [2]:
import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"]="2"
InĀ [3]:
# Properties
E = 50
nu = 0.3
a = 10
b = 10
f = 1
h = 1
1.3.2. Define PDE with Boundary ConditionsĀ¶
InĀ [4]:
def boundary_x(X, on_boundary):
return on_boundary and np.isclose(X[0], a)
def boundary_y(X, on_boundary):
return on_boundary and np.isclose(X[1], b)
def boundary_y0(X, on_boundary):
return on_boundary and np.isclose(X[1], 0)
def boundary_x0(X, on_boundary):
return on_boundary and np.isclose(X[0], 0)
InĀ [5]:
def pde(X, Y):
du_xx = dde.grad.hessian(Y, X, i = 0, j = 0, component = 0)
du_yy = dde.grad.hessian(Y, X, i = 1, j = 1, component = 0)
du_xy = dde.grad.hessian(Y, X, i = 0, j = 1, component = 0)
dv_xx = dde.grad.hessian(Y, X, i = 0, j = 0, component = 1)
dv_yy = dde.grad.hessian(Y, X, i = 1, j = 1, component = 1)
dv_xy = dde.grad.hessian(Y, X, i = 0, j = 1, component = 1)
force_eq_x = (du_xx + 0.5*(1 - nu)*du_yy + 0.5*(1 + nu)*dv_xy)*E/(1 - nu**2)
force_eq_y = (dv_yy + 0.5*(1 - nu)*dv_xx + 0.5*(1 + nu)*du_xy)*E/(1 - nu**2)
return [force_eq_x, force_eq_y]
InĀ [6]:
def BC_xy(X, Y):
du_x = dde.grad.jacobian(Y, X, i = 0, j = 0)
du_y = dde.grad.jacobian(Y, X, i = 0, j = 1)
dv_x = dde.grad.jacobian(Y, X, i = 1, j = 0)
dv_y = dde.grad.jacobian(Y, X, i = 1, j = 1)
sig_xx = (du_x + nu*dv_y)*E/(1 - nu**2)
sig_yy = (dv_y + nu*du_x)*E/(1 - nu**2)
sig_xy = (dv_x + du_y)*E/(1 + nu)/2
sig_ex = f * h * torch.cos(3.1415/(2*b)*X[:,1]).reshape(-1,1)
return sig_xx-sig_ex, sig_xy, sig_yy
1.3.3. Define Geometry and Implement Boundary ConditionĀ¶
InĀ [7]:
geom = dde.geometry.Rectangle(xmin=[0, 0], xmax=[a, b])
bc_Y0 = dde.DirichletBC(geom, lambda x:0, boundary_y0, component = 1)
bc_X0 = dde.DirichletBC(geom, lambda x:0, boundary_x0, component = 0)
bc_free_yy = dde.OperatorBC(geom, lambda x, y, _ : BC_xy(x,y)[2], boundary_y)
bc_free_yx = dde.OperatorBC(geom, lambda x, y, _ : BC_xy(x,y)[1], boundary_y)
bc_traction = dde.OperatorBC(geom, lambda x, y, _ : BC_xy(x,y)[0], boundary_x)
bc_free_xy = dde.OperatorBC(geom, lambda x, y, _ : BC_xy(x,y)[1], boundary_x)
InĀ [8]:
data = dde.data.PDE(geom,
pde,
[bc_traction, bc_free_xy, bc_free_yx, bc_free_yy, bc_Y0, bc_X0],
num_domain = 10000,
num_boundary = 1000,
num_test = 100,
train_distribution = 'LHS')
InĀ [9]:
plt.figure(figsize = (5,5))
plt.scatter(data.train_x_all[:,0], data.train_x_all[:,1], s = 0.1)
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.show()
1.3.4. Define Network and Hyper-parametersĀ¶
InĀ [10]:
layer_size = [2] + [64] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3)
1.3.5. Train (Adam Optimizer)Ā¶
InĀ [11]:
losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
1.3.6. Plot Displacement (Adam Optimzer)Ā¶
InĀ [12]:
color_legend = [[0, 0.182], [-0.06, 0.011], [-0.0022,1.0], [-0.15, 0.45]]
title = ['x-displacement ($u$)', 'y-displacement ($v$)', '$\sigma_{xx}$', '$\sigma_{yy}$']
InĀ [13]:
samples = geom.uniform_points(10201)
pde_disp = model.predict(samples)
plt.figure(figsize = (10, 4))
for idx in range(2):
plt.subplot(1,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = pde_disp[:, idx],
cmap = 'rainbow',
s = 5)
plt.clim(color_legend[idx])
plt.title(title[idx])
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
1.3.7. Train More (L-BFGS Optimizer)Ā¶
InĀ [14]:
dde.optimizers.config.set_LBFGS_options(maxiter = 5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
1.3.8. Plot Results (Adam + L-BFGS)Ā¶
InĀ [15]:
pde_disp = model.predict(samples)
plt.figure(figsize = (10, 4))
for idx in range(2):
plt.subplot(1,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = pde_disp[:, idx],
cmap = 'rainbow',
s = 5)
plt.clim(color_legend[idx])
plt.title(title[idx], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.savefig('PDE_disp.png', dpi=300, transparent = True)
plt.show()
1.3.9 Plot StressĀ¶
- Since the neural network's output is only x and y displacement, we should autograd to calculate stress
- If we assume that the Cauchy stress tensor components are linearly related to the von KƔrmƔn strains by Hooke's law, stress can be represented as:
$$
\sigma_{xx} = {E \over 1 - \nu^2}\left({\partial u \over \partial x} + \nu {\partial v \over \partial y} \right)\\\\
\sigma_{yy} = {E \over 1 - \nu^2}\left({\partial v \over \partial y} + \nu {\partial u \over \partial x} \right)
$$
InĀ [16]:
def check_stress(net, X):
X = torch.tensor(X, requires_grad = True)
disp = net(X)
u_pred = disp[:,0].reshape(-1,1)
v_pred = disp[:,1].reshape(-1,1)
du_xy = torch.autograd.grad(u_pred, X, torch.ones_like(u_pred),
retain_graph=True, create_graph=True, allow_unused=True)[0].detach().cpu().numpy()
dv_xy = torch.autograd.grad(v_pred, X, torch.ones_like(v_pred),
retain_graph=True, create_graph=True, allow_unused=True)[0].detach().cpu().numpy()
du_x, du_y = du_xy[:,0], du_xy[:,1]
dv_x, dv_y = dv_xy[:,0], dv_xy[:,1]
sig_xx = ((du_x + nu*dv_y)*E/(1 - nu**2)).reshape(-1,1)
sig_yy = ((dv_y + nu*du_x)*E/(1 - nu**2)).reshape(-1,1)
return sig_xx, sig_yy
InĀ [17]:
sig_x, sig_y = check_stress(model.net, samples)
pde_sig = np.hstack([sig_x, sig_y])
plt.figure(figsize = (10, 4))
for idx in range(2):
plt.subplot(1,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = pde_sig[:, idx],
cmap = 'rainbow',
s = 5)
plt.clim(color_legend[idx+2])
plt.title(title[idx+2], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.savefig('PDE_stress.png', dpi=300, transparent = True)
plt.show()
InĀ [20]:
diagonal_index = []
for i in range(u.shape[0]):
if loc[i,0] + loc[i,1] == 10:
diagonal_index.append(i)
diagonal_index = np.array(diagonal_index)
InĀ [21]:
diag_result_true = [u[diagonal_index], v[diagonal_index], stress[diagonal_index, 0], stress[diagonal_index, 1]]
diag_result_pde = [pde_disp[diagonal_index, 0], pde_disp[diagonal_index, 1], pde_sig[diagonal_index, 0], pde_sig[diagonal_index, 1]]
plt.figure(figsize = (9, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.plot(loc[diagonal_index,0], diag_result_true[idx], label = 'FEM')
plt.plot(loc[diagonal_index,0], diag_result_pde[idx], label = 'PDE')
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel(title[idx], fontsize = 14)
plt.xlim((0, a))
plt.legend(fontsize = 12)
plt.tight_layout()
plt.savefig('PDE_diag.png', dpi=300, transparent = True)
plt.show()
InĀ [23]:
Plate_data = pd.read_csv('./data/Plate_data.csv')
loc = Plate_data.iloc[:, 0:2].to_numpy()
u = Plate_data.iloc[:, 2:3].to_numpy()
v = Plate_data.iloc[:, 3:4].to_numpy()
stress = Plate_data.iloc[:, 4:7].to_numpy()
InĀ [24]:
random.seed(1)
tr_sample = random.sample(range(0,Plate_data.shape[0]),5000)
x_sample = loc[tr_sample,:]
u_sample = u[tr_sample,:]
v_sample = v[tr_sample,:]
InĀ [25]:
observe_u = dde.icbc.PointSetBC(x_sample, u_sample, component=0)
observe_v = dde.icbc.PointSetBC(x_sample, v_sample, component=1)
2.2. Define GeometryĀ¶
InĀ [26]:
geom = dde.geometry.Rectangle(xmin=[0, 0], xmax=[a, b])
data = dde.data.PDE(geom,
None,
[observe_u, observe_v],
num_domain = 0,
num_boundary = 0,
num_test = 100)
InĀ [27]:
plt.figure(figsize = (5,5))
plt.scatter(data.train_x_bc[:,0], data.train_x_bc[:,1], s = 1)
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.show()
2.3. Define Network and Hyper-parametersĀ¶
InĀ [28]:
layer_size = [2] + [64] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3)
2.4. Train (Adam Optimizer)Ā¶
InĀ [29]:
losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = False)
2.5. Train More (L-BFGS Optimizer)Ā¶
InĀ [30]:
dde.optimizers.config.set_LBFGS_options(maxiter = 5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
2.6. Plot Results (Adam + L-BFGS)Ā¶
InĀ [31]:
big_disp = model.predict(samples)
sig_x, sig_y = np.array(check_stress(model.net, samples))
big_sig = np.hstack([sig_x, sig_y])
big_result = np.hstack([big_disp, big_sig])
plt.figure(figsize = (10, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = big_result[:,idx],
cmap = 'rainbow',
s = 5)
plt.title(title[idx], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.savefig('Bigdata_disp_stress.png', dpi=300, transparent = True)
plt.show()
InĀ [33]:
diag_result_true = [u[diagonal_index], v[diagonal_index], stress[diagonal_index, 0], stress[diagonal_index, 1]]
diag_result_big = [big_disp[diagonal_index, 0], big_disp[diagonal_index, 1], big_sig[diagonal_index, 0], big_sig[diagonal_index, 1]]
plt.figure(figsize = (9, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.plot(loc[diagonal_index,0], diag_result_true[idx], label = 'FEM')
plt.plot(loc[diagonal_index,0], diag_result_big[idx], label = 'Big data')
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel(title[idx], fontsize = 14)
plt.xlim((0, a))
plt.legend(fontsize = 12)
plt.tight_layout()
plt.savefig('Big_diag.png', dpi=300, transparent = True)
plt.show()
InĀ [35]:
random.seed(1)
tr_sample = random.sample(range(0,Plate_data.shape[0]),100)
x_sample = loc[tr_sample,:]
u_sample = u[tr_sample,:]
v_sample = v[tr_sample,:]
InĀ [36]:
observe_u = dde.icbc.PointSetBC(x_sample, u_sample, component=0)
observe_v = dde.icbc.PointSetBC(x_sample, v_sample, component=1)
3.2. Define GeometryĀ¶
InĀ [37]:
geom = dde.geometry.Rectangle(xmin=[0, 0], xmax=[a, b])
data = dde.data.PDE(geom,
None,
[observe_u, observe_v],
num_domain = 0,
num_boundary = 0,
num_test = 100)
InĀ [38]:
plt.figure(figsize = (5,5))
plt.scatter(data.train_x_bc[:,0], data.train_x_bc[:,1], s = 1)
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.show()
3.3. Define Network and Hyper-parametersĀ¶
InĀ [39]:
layer_size = [2] + [64] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3)
3.4. Train (Adam Optimizer)Ā¶
InĀ [40]:
losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = False)
3.5. Train More (L-BFGS Optimizer)Ā¶
InĀ [41]:
dde.optimizers.config.set_LBFGS_options(maxiter = 5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
3.6. Plot Results (Adam + L-BFGS)Ā¶
InĀ [42]:
small_disp = model.predict(samples)
sig_x, sig_y = np.array(check_stress(model.net, samples))
small_sig = np.hstack([sig_x, sig_y])
small_result = np.hstack([small_disp, small_sig])
plt.figure(figsize = (10, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = small_result[:,idx],
cmap = 'rainbow',
s = 5)
plt.title(title[idx], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.savefig('Smalldata_disp_stress.png', dpi=300, transparent = True)
plt.show()
InĀ [44]:
diag_result_true = [u[diagonal_index], v[diagonal_index], stress[diagonal_index, 0], stress[diagonal_index, 1]]
diag_result_small = [small_disp[diagonal_index, 0], small_disp[diagonal_index, 1], small_sig[diagonal_index, 0], small_sig[diagonal_index, 1]]
plt.figure(figsize = (9, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.plot(loc[diagonal_index,0], diag_result_true[idx], label = 'FEM')
plt.plot(loc[diagonal_index,0], diag_result_small[idx], label = 'Small data')
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel(title[idx], fontsize = 14)
plt.xlim((0, a))
plt.legend(fontsize = 12)
plt.tight_layout()
plt.savefig('Small_diag.png', dpi=300, transparent = True)
plt.show()
4. PINN with Small DataĀ¶
InĀ [46]:
data = dde.data.PDE(geom,
pde,
[bc_traction, bc_free_xy, bc_free_yx, bc_free_yy, bc_Y0, bc_X0, observe_u, observe_v],
num_domain = 100,
num_boundary = 100,
num_test = 100,
train_distribution = 'LHS')
InĀ [47]:
plt.figure(figsize = (5,5))
plt.scatter(data.train_x_all[:,0], data.train_x_all[:,1], s = 1)
plt.xlabel('x (mm)')
plt.ylabel('y (mm)')
plt.show()
4.1. Define Network and Hyper-parametersĀ¶
InĀ [48]:
layer_size = [2] + [64] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3, loss_weights=[1, 1, 1, 1, 1, 1, 1, 1, 9, 9])
4.2 Train (Adam Optimizer)Ā¶
InĀ [49]:
losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = False)
4.3. Train More (L-BFGS Optimizer)Ā¶
InĀ [50]:
dde.optimizers.config.set_LBFGS_options(maxiter = 5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
4.4. Plot Results (Adam + L-BFGS)Ā¶
InĀ [51]:
pde_small_disp = model.predict(samples)
sig_x, sig_y = np.array(check_stress(model.net, samples))
pde_small_sig = np.hstack([sig_x, sig_y])
pde_small_result = np.hstack([pde_small_disp, pde_small_sig])
plt.figure(figsize = (10, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.scatter(samples[:, 0],
samples[:, 1],
c = pde_small_result[:,idx],
cmap = 'rainbow',
s = 5)
plt.title(title[idx], fontsize = 15)
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel('y (mm)', fontsize = 14)
plt.axis('square')
plt.xlim((0, a))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.savefig('PDEdata_disp_stress.png', dpi=300, transparent = True)
plt.show()
InĀ [53]:
diag_result_true = [u[diagonal_index], v[diagonal_index], stress[diagonal_index, 0], stress[diagonal_index, 1]]
diag_result_pde_small = [pde_small_disp[diagonal_index, 0], pde_small_disp[diagonal_index, 1], pde_small_sig[diagonal_index, 0], pde_small_sig[diagonal_index, 1]]
plt.figure(figsize = (9, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.plot(loc[diagonal_index,0], diag_result_true[idx], label = 'FEM')
plt.plot(loc[diagonal_index,0], diag_result_pde_small[idx], label = 'PDE+small data')
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel(title[idx], fontsize = 14)
plt.xlim((0, a))
plt.legend(fontsize = 12)
plt.tight_layout()
plt.savefig('PDE+small_diag.png', dpi=300, transparent = True)
plt.show()
4.5. Plot with All ResultsĀ¶
InĀ [64]:
plt.figure(figsize = (9, 8))
for idx in range(4):
plt.subplot(2,2,idx+1)
plt.plot(loc[diagonal_index,0], diag_result_true[idx],linewidth = 2, c = 'k' ,label = 'FEM')
plt.plot(loc[diagonal_index,0], diag_result_pde[idx], '--' ,linewidth = 2 ,label = 'PDE')
plt.plot(loc[diagonal_index,0], diag_result_big[idx], '-.' ,linewidth = 2, label = 'Big data')
plt.plot(loc[diagonal_index,0], diag_result_small[idx], ':' ,linewidth = 2 , label = 'Small data')
plt.plot(loc[diagonal_index,0], diag_result_pde_small[idx], '.',linewidth = 2 , label = 'PDE+small data')
plt.xlabel('x (mm)', fontsize = 14)
plt.ylabel(title[idx], fontsize = 14)
plt.xlim((0, a))
plt.legend(fontsize = 11)
plt.tight_layout()
plt.savefig('Total_diag.png', dpi=300, transparent = True)
plt.show()
InĀ [3]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')